// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtPVLDFileFormat.C                            //
// ************************************************************************* //

#include <avtPVLDFileFormat.h>

#include <string>
#include <sstream>
#include <algorithm>
#include <stdexcept>

#include <vtkCommand.h>
#include <vtkIntArray.h>
#include <vtkFloatArray.h>
#include <vtkCellType.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>
#include <avtMaterial.h>
#include <avtVariableCache.h>

#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>
#include <DebugStream.h>

#include <dbginfo.h>
#include <string>
#include <sstream>
#include <algorithm>
#include <stdexcept>

using     std::vector;
using     std::string;

//#define DEBUG_PRINT

static PVLD_Reader *
InstantiateReader()
{
    PVLD_Reader *r = NULL;
    if(getenv("ELEMENT_PARTITION_VLD") == NULL)
        r = new PVLD_Part_Reader;
    else
        r = new PVLD_Reader;

    return r;
}


// ****************************************************************************
//  Method: avtPVLDFileFormat constructor
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
//  Modifications:
//   Brad Whitlock, Mon Jun  2 18:40:34 PDT 2014
//   Use a helper function to instantiate the preader.
//
// ****************************************************************************

avtPVLDFileFormat::avtPVLDFileFormat(const char *filename)
    : avtSTMDFileFormat(&filename, 1)
{
    debug5 << "creating reader( \"" << filename << "\")...\n";

    {
        const char *spt = getenv("ADD_MISSING_PARTS_VLD");
        add_missing_parts_ = spt!=NULL;
    }

    enablePVLD = (getenv("PVLDDISABLENEW") == NULL);
    hasTOCread_ = false;
    preader_ = InstantiateReader();
    preader_->SetFileName( filename );
}

avtPVLDFileFormat::~avtPVLDFileFormat()
{
    delete preader_;
}

// ****************************************************************************
//  Method: avtBOVFileFormat::ActivateTimestep
//
//  Purpose:
//      Activates the current timestep, prompting the header to be read.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Mon Sep 17 22:07:29 PST 2012
//
// ****************************************************************************

void
avtPVLDFileFormat::ActivateTimestep(void)
{
    debug1 << "calling avtPVLDFileFormat::ActivateTimestep()...\n";
    if( !hasTOCread_ )
    {
        OutputMemoryInfo( "Before ReadTOC()" );
        ReadTOC();
        hasTOCread_=true;
        OutputMemoryInfo( "After ReadTOC()" );
    }
}


// ****************************************************************************
//  Method: avtPVLDFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
// Modifications:
//   Brad Whitlock, Mon Jun  2 18:40:34 PDT 2014
//   Totally delete the preader when we free up resources so all vectors it
//   contains will be deleted automatically.
//
// ****************************************************************************

void
avtPVLDFileFormat::FreeUpResources(void)
{
    debug1 << "calling avtPVLDFileFormat::FreeUpResources()...\n";
    OutputMemoryInfo( "Before FreeUpResources()" );

    delete preader_;
    hasTOCread_ = false;
    preader_ = InstantiateReader();
    preader_->SetFileName( filenames[0] );

    OutputMemoryInfo( "After FreeUpResources()" );
}

// ****************************************************************************
//  Method: avtPVLDFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
//  Modifications:
//    Brad Whitlock, Mon Jun  2 18:42:20 PDT 2014
//    Store the PVLD title into the database comment.
//
//    Brad Whitlock, Wed Jun 18 17:04:54 PDT 2014
//    Read the TOC if it has not been read yet.
//
// ****************************************************************************

void
avtPVLDFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    debug1 << "calling avtPVLDFileFormat::PopulateDatabaseMetaData()...\n";
    TRY
    {
        // Read the TOC in case it has not been read and ActivateTimestep was not
        // called, as it seems can be the case on the mdserver.
        if( !hasTOCread_ )
        {
            OutputMemoryInfo( "Before ReadTOC()" );
            ReadTOC();
            hasTOCread_=true;
            OutputMemoryInfo( "After ReadTOC()" );
        }

        md->SetDatabaseComment(preader_->GetTitle());

        AddSolidMesh( md );
        AddShellMesh( md );
        AddBeamMesh(  md );
        AddSurfaceMesh( md );
        AddSphMesh( md );
        AddTiedSetMesh( md );
        AddContactMesh( md );

        AddSolidVariables( md );
        AddShellVariables( md );
        AddBeamVariables(  md );
        AddSurfaceVariables( md );
        AddSphVariables( md );
        AddTiedSetVariables( md );
        AddContactVariables( md );

        AddSolidMaterial( md );
        AddShellMaterial( md );
        AddBeamMaterial(  md );
        AddSurfaceMaterial( md );
        AddSphMaterial( md );

        AddSolidIndexVariables( md );
    }
    CATCH2( std::exception, e )
    {
        string msg = e.what();
        msg += "Failed in avtPVLDFileFormat::PopulateDatabaseMetaData()\n";
        EXCEPTION1(InvalidDBTypeException, e.what() );
    }
    ENDTRY
}

// ****************************************************************************
// Method: avtPVLDFileFormat::CreateGlobalNodeIds
//
// Purpose:
//   Create global node ids for the domain.
//
// Arguments:
//   domain   : The domain
//   meshname : The name of the mesh
//   vmap     : A vector containing the global node ids.
//
// Programmer: Brad Whitlock
// Creation:   Mon Jun  2 18:23:54 PDT 2014
//
// Modifications:
//
// ****************************************************************************

void
avtPVLDFileFormat::CreateGlobalNodeIds(int domain, const char *meshname, const vector<int> &vmap)
{
    if(enablePVLD)
    {
        // Use the vmap from the reader to expose global node ids that VisIt can use
        // for ghost zone creation.

// NOTE: This assumes that there were no missing material mesh cells added by
//       AppendMissingMaterialMesh.

        vtkIntArray *arr = vtkIntArray::New();
        arr->SetNumberOfTuples(vmap.size());
        arr->SetName("avtGlobalNodeId");
        memcpy((int *)arr->GetVoidPointer(0), &vmap[0], vmap.size() * sizeof(int));

        //
        // Cache this VTK object but in the VoidRefCache, not the VTK cache
        // so that it can be obtained through the GetAuxiliaryData call
        //
        debug1 << "Storing global node ids for " << meshname << " domain " << domain << endl;
        void_ref_ptr vr = void_ref_ptr(arr, avtVariableCache::DestructVTKObject);
        cache->CacheVoidRef(meshname, AUXILIARY_DATA_GLOBAL_NODE_IDS, timestep,
                            domain, vr);
    }
}

// ****************************************************************************
//  Method: avtPVLDFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
//  Modifications:
//    Brad Whitlock, Thu May 29 16:22:23 PST 2014
//    I added global node ids for so VisIt can automatically generate ghost
//    cells and remove artifacts at domain boundaries.
//
// ****************************************************************************

vtkDataSet *
avtPVLDFileFormat::GetMesh(int domain, const char *meshname)
{
    debug1 << "calling avtPVLDFileFormat::GetMesh("
           << domain << ", " << meshname << ")\n";
    TRY
    {
        vector<int>   ele;
        vector<float> crd;
        if( PVLD_Reader::solid_name == meshname )
        {
            preader_->CheckNumberOfEngines( preader_->GetNumOfSolidBlocks() );
            preader_->ReadSolidBlockMesh( domain, crd, ele );
            preader_->AppendMissingMaterialMesh( PVLD_Reader::solid_elmt_type, domain, crd, ele );

            if( ele.size()>0 )
            {
                vtkDataSet *ds = GenerateSolidMesh( crd, ele );
                CreateGlobalNodeIds(domain, meshname, preader_->GetSolidBlockMeshMap(domain));
                return ds;
            }
        }
        else if( PVLD_Reader::shell_name == meshname )
        {
            preader_->CheckNumberOfEngines( preader_->GetNumOfShellBlocks() );
            preader_->ReadShellBlockMesh( domain, crd, ele );
            preader_->AppendMissingMaterialMesh( PVLD_Reader::shell_elmt_type, domain, crd, ele );

            if( ele.size()>0 )
            {
                vtkDataSet *ds = GenerateShellMesh( crd, ele );
                CreateGlobalNodeIds(domain, meshname, preader_->GetShellBlockMeshMap(domain));
                return ds;
            }
        }
        else if( PVLD_Reader::beam_name == meshname )
        {
            preader_->CheckNumberOfEngines( preader_->GetNumOfBeamBlocks() );
            preader_->ReadBeamBlockMesh( domain, crd, ele );
            preader_->AppendMissingMaterialMesh( PVLD_Reader::beam_elmt_type, domain, crd, ele );
            if( ele.size()>0 ) return GenerateBeamMesh( crd, ele );

            if( ele.size()>0 )
            {
                vtkDataSet *ds = GenerateBeamMesh( crd, ele );
                CreateGlobalNodeIds(domain, meshname, preader_->GetBeamBlockMeshMap(domain));
                return ds;
            }
        }
        else if( PVLD_Reader::surface_name == meshname )
        {
            preader_->CheckNumberOfEngines( preader_->GetNumOfSurfaceBlocks() );
            preader_->ReadSurfaceBlockMesh( domain, crd, ele );

            if( ele.size()>0 )
            {
                vtkDataSet *ds = GenerateShellMesh( crd, ele );
                CreateGlobalNodeIds(domain, meshname, preader_->GetSurfaceBlockMeshMap(domain));
                return ds;
            }
        }
        else if( PVLD_Reader::sph_name == meshname )
        {
            TRY
            {
                preader_->ReadSphBlockMesh( domain, crd, ele );
                preader_->AppendMissingMaterialMesh( PVLD_Reader::sph_elmt_type, domain, crd, ele );
                if( ele.size()>0 )
                    return GenerateSphMesh( crd, ele );
            }
            CATCH(std::exception)
            {
                char msg[512];
                snprintf(msg, 512, "SPH data may not be available at time step %d", timestep);
                EXCEPTION1(VisItException, msg);
            }
            ENDTRY
        }
        else {
            vector<string> vs;
            PVLD_Reader::DecomposeNames( meshname, vs );

            if( vs[0] == PVLD_Reader::tiedset_name )
            {
                if( vs[1] == "Segment" )
                {
                    preader_->ReadTiedSetBlockMesh( domain, crd, ele );
                    if( crd.size()>0 ) return GenerateShellMesh( crd, ele );
                }
                else if( vs[1] == "Master" )
                {
                    preader_->ReadTiedSetManagerBlockMesh( domain, crd );
                    if( crd.size()>0 ) return GeneratePointMesh( crd );
                }
                else if( vs[1] == "Slave" )
                {
                    preader_->ReadTiedSetWorkerBlockMesh( domain, crd );
                    if( crd.size()>0 ) return GeneratePointMesh( crd );
                }
                else if( vs[1] == "Bond" )
                {
                    vector<float> crd2;
                    preader_->ReadTiedSetManagerBlockMesh( domain, crd );
                    preader_->ReadTiedSetWorkerBlockMesh( domain, crd2 );
                    if( crd.size()>0 ) return GenerateBondMesh( crd, crd2 );
                }
            }
            else if( vs[0] == PVLD_Reader::contact_name )
            {
                if( vs[1] == "Segment" )
                {
                    preader_->ReadContactBlockMesh( domain, crd, ele );
                    if( crd.size()>0 ) return GenerateShellMesh( crd, ele );
                }
                else if( vs[1] == "Master" )
                {
                    preader_->ReadContactManagerBlockMesh( domain, crd );
                    if( crd.size()>0 ) return GeneratePointMesh( crd );
                }
                else if( vs[1] == "Slave" )
                {
                    preader_->ReadContactWorkerBlockMesh( domain, crd );
                    if( crd.size()>0 ) return GeneratePointMesh( crd );
                }
                else if( vs[1] == "Bond" )
                {
                    vector<float> crd2;
                    preader_->ReadContactManagerBlockMesh( domain, crd );
                    preader_->ReadContactWorkerBlockMesh( domain, crd2 );
                    if( crd.size()>0 ) return GenerateBondMesh( crd, crd2 );
                }
            }
        }
    }
    CATCH2( std::exception, e )
    {
        string msg = e.what();
        msg += "Failed in avtPVLDFileFormat::GetMesh()\n";
        EXCEPTION1(InvalidDBTypeException, e.what() );
    }
    ENDTRY
    return NULL;
}


// ****************************************************************************
//  Method: avtPVLDFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
// ****************************************************************************

vtkDataArray *
avtPVLDFileFormat::GetVar(int domain, const char *varname)
{
    debug1 << "calling avtPVLDFileFormat::GetVar("
           << domain << ", " << varname << ")\n";
    TRY
    {
        vector<string> vecname;
        PVLD_Reader::DecomposeNames( varname, vecname );
        string& mname = vecname[0];
        string& vname = vecname[1];

        vector<int> dims;
        vector<float> data;
        if( mname == PVLD_Reader::solid_name )
        {
            if( vecname.size()==2 )
                preader_->ReadSolidBlockData( vname.c_str(), domain, dims, data );
            else
                preader_->ReadSolidBlockIndexVariable( vname, vecname[2], domain, dims, data );
            preader_->AppendMissingMaterialData( PVLD_Reader::solid_elmt_type, domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
        else if( mname == PVLD_Reader::shell_name )
        {
            if( vecname.size()==2 )
                preader_->ReadShellBlockData( vname.c_str(), domain, dims, data );
            else
                preader_->ReadShellBlockIndexVariable( vname, vecname[2], domain, dims, data );
            preader_->AppendMissingMaterialData( PVLD_Reader::shell_elmt_type, domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
        else if( mname == PVLD_Reader::beam_name )
        {
            if( vecname.size()==2 )
                preader_->ReadBeamBlockData( vname.c_str(), domain, dims, data );
            else
                //preader_->ReadBeamBlockIndexVariable( vname, vecname[2], domain, dims, data );
                throw std::runtime_error( "No IndexVariable support in Beam!\n" );
            preader_->AppendMissingMaterialData( PVLD_Reader::beam_elmt_type, domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
        else if( mname == PVLD_Reader::surface_name )
        {
            preader_->ReadSurfaceBlockData( vname.c_str(), domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
        else if( mname == PVLD_Reader::sph_name )
        {
            TRY
            {
                preader_->ReadSphBlockData( vname.c_str(), domain, dims, data );
                preader_->AppendMissingMaterialData( PVLD_Reader::sph_elmt_type, domain, dims, data );
                if( data.size()>0 ) return GenerateVariable( dims, data );
            }
            CATCH(std::exception)
            {
                char msg[512];
                snprintf(msg, 512, "SPH data may not be available at time step %d", timestep);
                EXCEPTION1(VisItException, msg);
            }
            ENDTRY
        }
        else if( mname == PVLD_Reader::tiedset_name )
        {
            string& vn = vecname[2];
            preader_->ReadTiedSetBlockData( vn.c_str(), domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
        else if( mname == PVLD_Reader::contact_name )
        {
            string& vn = vecname[2];
            preader_->ReadContactBlockData( vn.c_str(), domain, dims, data );
            if( data.size()>0 ) return GenerateVariable( dims, data );
        }
    }
    CATCH2( std::exception, e )
    {
        string msg = e.what();
        msg += "Failed in avtPVLDFileFormat::GetVar()\n";
        EXCEPTION1(InvalidVariableException, e.what() );
    }
    ENDTRY
    return NULL;
}


// ****************************************************************************
//  Method: avtPVLDFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: hpan -- generated by xml2avt
//  Creation:   Tue Sep 25 14:55:38 PST 2012
//
// ****************************************************************************

vtkDataArray *
avtPVLDFileFormat::GetVectorVar(int domain, const char *varname)
{
    debug1 << "calling avtPVLDFileFormat::GetVectorVar( "
           << domain << ", " << varname << ")\n";
    TRY
    {
        string mname, vname;
        PVLD_Reader::DecomposeNames( varname, mname, vname );

        vector<int> dims;
        vector<float> data;
        if( mname == PVLD_Reader::solid_name )
        {
            preader_->ReadSolidBlockData( vname.c_str(), domain, dims, data );
            preader_->AppendMissingMaterialData( PVLD_Reader::solid_elmt_type, domain, dims, data );
            if( data.size()>0 )
            {
                if( dims.size()==2 && dims[1]==6 )
                    return GenerateDyna3dStressTensor( dims, data );
                else
                    return GenerateVariable( dims, data );
            }
        }
        else if( mname == PVLD_Reader::shell_name )
        {
            preader_->ReadShellBlockData( vname.c_str(), domain, dims, data );
            preader_->AppendMissingMaterialData( PVLD_Reader::shell_elmt_type, domain, dims, data );
            if( data.size()>0 )
            {
                if( dims.size()==2 && dims[1]==6 )
                    return GenerateDyna3dStressTensor( dims, data );
                else
                    return GenerateVariable( dims, data );
            }
        }
        else if( mname == PVLD_Reader::beam_name )
        {
            preader_->ReadBeamBlockData( vname.c_str(), domain, dims, data );
            preader_->AppendMissingMaterialData( PVLD_Reader::beam_elmt_type, domain, dims, data );
            if( data.size()>0 )
            {
                if( dims.size()==2 && dims[1]==6 )
                    return GenerateDyna3dStressTensor( dims, data );
                else
                    return GenerateVariable( dims, data );
            }
        }
        else if( mname == PVLD_Reader::surface_name )
        {
            preader_->ReadSurfaceBlockData( vname.c_str(), domain, dims, data );
            if( data.size()>0 )
            {
                if( dims.size()==2 && dims[1]==6 )
                    return GenerateDyna3dStressTensor( dims, data );
                else
                    return GenerateVariable( dims, data );
            }
        }
        else if( mname == PVLD_Reader::sph_name )
        {
            TRY
            {
                preader_->ReadSphBlockData( vname.c_str(), domain, dims, data );
                preader_->AppendMissingMaterialData( PVLD_Reader::sph_elmt_type, domain, dims, data );
                if( data.size()>0 )
                {
                    if( dims.size()==2 && dims[1]==6 )
                        return GenerateSphStressTensor( dims, data );
                    else
                        return GenerateVariable( dims, data );
                }
            }
            CATCH(std::exception)
            {
                char msg[512];
                snprintf(msg, 512, "SPH data may not be available at time step %d", timestep);
                EXCEPTION1(VisItException, msg);
            }
            ENDTRY
        }
    }
    CATCH2( std::exception, e )
    {
        string msg = e.what();
        msg += "Failed in avtPVLDFileFormat::GetVectorVar()\n";
        EXCEPTION1(InvalidDBTypeException, e.what() );
    }
    ENDTRY
    return NULL;
}




void*
avtPVLDFileFormat::
GetAuxiliaryData( const char *varname, int domain,
                  const char *type, void *args,
                  DestructorFunction &df )
{
    debug1 << "calling avtPVLDFileFormat::GetAuxiliaryData( "
           << varname << ", " << domain << ", " << type << ", ...)\n";
    void* retval=NULL;

    TRY
    {
        if(strcmp(type, AUXILIARY_DATA_MATERIAL) == 0)
        {
            debug1 << "Retrieve AUXILIARY_DATA_MATERIAL\n";
            string mname, vname;
            PVLD_Reader::DecomposeNames( varname, mname, vname, "_" );

            vector<int> mat;
            vector<int> ddms;
            vector<float> data;
            if( mname == PVLD_Reader::solid_name )
            {
                preader_->ReadSolidMaterial();
                mat = preader_->GetSolidMaterial();
                const vector<int> &msm = preader_->GetMissingSolidMaterial();
                mat.insert( mat.end(), msm.begin(), msm.end() );
                //preader_->ReadSolidBlockData( "Material", domain, ddms, data );
                preader_->ReadSolidBlockMaterial( domain, ddms, data );
                preader_->AppendMissingMaterialMaterial( PVLD_Reader::solid_elmt_type, domain, ddms, data );
            }
            else if( mname == PVLD_Reader::shell_name )
            {
                preader_->ReadShellMaterial();
                mat = preader_->GetShellMaterial();
                const vector<int> &msm = preader_->GetMissingShellMaterial();
                mat.insert( mat.end(), msm.begin(), msm.end() );
                //preader_->ReadShellBlockData( "Material", domain, ddms, data );
                preader_->ReadShellBlockMaterial( domain, ddms, data );
                preader_->AppendMissingMaterialMaterial( PVLD_Reader::shell_elmt_type, domain, ddms, data );
            }
            else if( mname == PVLD_Reader::beam_name )
            {
                preader_->ReadBeamMaterial();
                mat = preader_->GetBeamMaterial();
                const vector<int> &msm = preader_->GetMissingBeamMaterial();
                mat.insert( mat.end(), msm.begin(), msm.end() );
                preader_->ReadBeamBlockMaterial( domain, ddms, data );
                preader_->AppendMissingMaterialMaterial( PVLD_Reader::beam_elmt_type, domain, ddms, data );
            }
            else if( mname == PVLD_Reader::surface_name )
            {
                preader_->ReadSurfaceMaterial();
                mat = preader_->GetSurfaceMaterial();
                //preader_->ReadSurfaceBlockData( "Material", domain, ddms, data );
                preader_->ReadSurfaceBlockMaterial( domain, ddms, data );
            }
            else if( mname == PVLD_Reader::sph_name )
            {
                TRY
                {
                    preader_->ReadSphMaterial(false);
                    mat = preader_->GetSphMaterial();
                    const vector<int> &msm = preader_->GetMissingSphMaterial();
                    mat.insert( mat.end(), msm.begin(), msm.end() );
                    //preader_->ReadSphBlockData( "Material", domain, ddms, data );
                    preader_->ReadSphBlockMaterial( domain, ddms, data );
                    preader_->AppendMissingMaterialMaterial( PVLD_Reader::sph_elmt_type, domain, ddms, data );
                }
                CATCH(std::exception)
                {
                    char msg[512];
                    snprintf(msg, 512, "SPH material data may not be available at time step %d", timestep);
                    EXCEPTION1(VisItException, msg);
                }
                ENDTRY
            }

            if( data.size()>0  )
            {
                if( mat.size()>0 )
                {
                    char** mnames= new char*[ mat.size() ];
                    int cnt=0;

                    for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                    {
                        int m = *iter;
                        const string& mt = preader_->GetMaterialName( m );
                        mnames[cnt] = new char[ mt.length()+1 ];
                        strcpy( mnames[cnt], mt.c_str() );
                        ++cnt;
                    }
                    vector<int> mid( data.size() );
                    for( size_t i=0; i<data.size(); i++ )
                        mid[i] = int( data[i] );

                    int ndims=1;
                    int dims[1];
                    dims[0]=(int)mid.size();
                    avtMaterial *amat = new avtMaterial(
                        (int)mat.size(), &mat[0],
                        mnames,
                        ndims, dims,
                        0, &mid[0],
                        0, // length of mix arrays
                        0, // mix_mat array
                        0, // mix_next array
                        0, // mix_zone array
                        0  // mix_vf array
                    );
                    for( size_t i=0; i<mat.size(); i++ )
                        delete [] mnames[i];
                    delete [] mnames;

                    retval = amat;
                    df = avtMaterial::Destruct;
                }
            }
        }
    }
    CATCH2( std::exception, e )
    {
        string msg = e.what();
        msg += "Failed in avtPVLDFileFormat::GetAuxiliaryData()\n";
        EXCEPTION1(InvalidDBTypeException, msg.c_str() );
    }
    ENDTRY
    return retval;
}




void
avtPVLDFileFormat::ReadTOC()
{
    TRY
    {
        preader_->ReadTOC();
        preader_->ReadHistoryTOC();
        debug1 << "add_missing_parts_ = " << add_missing_parts_ << "\n";
        preader_->ReadMaterialType( add_missing_parts_ );
    }
    CATCH2( std::exception, e )
    {
        debug1 << "Exception caught in avtVelodyneFileFormat::ReadTOC(): \n\t"
               << e.what() << ".\n";
        EXCEPTION1(InvalidDBTypeException, e.what() );
    }
    ENDTRY
#ifdef DEBUG_PRINT
    debug5 << "ncycle=" << preader_->GetCycle()
           << "\t time=" << preader_->GetTime()
           << "\t nmmat=" << preader_->GetNumOfMaterials()
           << "\n";

    debug5 << "HasNode=" << preader_->HasNode() << "\n";
    if( preader_->HasNode() )
    {
        debug5 << "\t NumNodes=" << preader_->GetNumOfNodes() << "\n"
               << "\t nNodeVariables=" << preader_->GetNumOfNodeVariables() << "\n";
        for( int i=0; i<preader_->GetNumOfNodeVariables(); i++ )
        {
            const string& name = preader_->GetNodeVariableName(i);
            const string& dims = preader_->GetNodeVariableDims(i);
            debug5 << "\t data_" << i << ": " << name << "  [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
    }

    debug5 << "HasSolid=" << preader_->HasSolid() << "\n";
    if( preader_->HasSolid() )
    {
        debug5 << "\t nSolidBlocks=" << preader_->GetNumOfSolidBlocks() << "\n"
               << "\t nSolidVariables=" << preader_->GetNumOfSolidVariables() << "\n";
        for( int i=0; i<preader_->GetNumOfSolidVariables(); i++ )
        {
            const string& name = preader_->GetSolidVariableName(i);
            const vector<int>& dims = preader_->GetSolidVariableDims(i);
            debug5 << "\t data_" << i << ": " << name << "  [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
        if( preader_->HasSolidHistoryVariables() )
        {
            debug5 << "\t MaxSolidHistoryVariables=" <<  preader_->MaxSolidHistoryVariables() << "\n";
            const vector<int>& dims = preader_->GetSolidHistoryPartitions();
            debug5 << "\t hvpart = [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
    }

    debug5 << "HasShell=" << preader_->HasShell() << "\n";
    if( preader_->HasShell() )
    {
        debug5 << "\t nShellBlocks=" << preader_->GetNumOfShellBlocks() << "\n"
               << "\t nShellVariables=" << preader_->GetNumOfShellVariables() << "\n";
        for( int i=0; i<preader_->GetNumOfShellVariables(); i++ )
        {
            const string& name = preader_->GetShellVariableName(i);
            const vector<int>& dims = preader_->GetShellVariableDims(i);
            debug5 << "\t data_" << i << ": " << name << "  [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
        if( preader_->HasShellHistoryVariables() )
        {
            debug5 << "\t MaxShellHistoryVariables=" <<  preader_->MaxShellHistoryVariables() << "\n";
            const vector<int>& dims = preader_->GetShellHistoryPartitions();
            debug5 << "\t hvpart = [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
    }

    debug5 << "HasSurface=" << preader_->HasSurface() << "\n";
    if( preader_->HasSurface() )
    {
        debug5 << "\t nSurfaceBlocks=" << preader_->GetNumOfSurfaceBlocks() << "\n"
               << "\t nSurfaceVariables=" << preader_->GetNumOfSurfaceVariables() << "\n";
        for( int i=0; i<preader_->GetNumOfSurfaceVariables(); i++ )
        {
            const string& name = preader_->GetSurfaceVariableName(i);
            const vector<int>& dims = preader_->GetSurfaceVariableDims(i);
            debug5 << "\t data_" << i << ": " << name << "  [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
    }

    debug5 << "HasSPH=" << preader_->HasSPH() << "\n";
    if( preader_->HasSPH() )
    {
        debug5 << "\t nSphBlocks=" << preader_->GetNumOfSphBlocks() << "\n"
               << "\t nSphVariables=" << preader_->GetNumOfSphVariables() << "\n";
        for( int i=0; i<preader_->GetNumOfSphVariables(); i++ )
        {
            const string& name = preader_->GetSphVariableName(i);
            const vector<int>& dims = preader_->GetSphVariableDims(i);
            debug5 << "\t data_" << i << ": " << name << "  [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
        if( preader_->HasSphHistoryVariables() )
        {
            debug5 << "\t MaxSphHistoryVariables=" <<  preader_->MaxSphHistoryVariables() << "\n";
            const vector<int>& dims = preader_->GetSphHistoryPartitions();
            debug5 << "\t hvpart = [";
            std::for_each( dims.begin(), dims.end(),
                           [&]( const int n )
            {
                debug5 << "  " << n;
            } );
            debug5 << "  ]\n";
        }
    }

    debug5 << "HasTiedSet=" << preader_->HasTiedSet() << "\n";
    if( preader_->HasTiedSet() )
    {
        debug5 << "\t nTiedSetBlocks=" << preader_->GetNumOfTiedSetBlocks() << "\n"
               << "\t nTiedSetVariables=" << preader_->GetNumOfTiedSetVariables() << "\n";
    }

    debug5 << "HasContact=" << preader_->HasContact() << "\n";
    if( preader_->HasContact() )
    {
        debug5 << "\t nContactBlocks=" << preader_->GetNumOfContactBlocks() << "\n"
               << "\t nContactVariables=" << preader_->GetNumOfContactVariables() << "\n";
    }
#endif
}


void
avtPVLDFileFormat::AddSolidMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasSolid() )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::solid_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 3;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfSolidBlocks();
        md->Add(mmd);
    }
}

void
avtPVLDFileFormat::AddShellMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasShell() )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::shell_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 2;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfShellBlocks();
        md->Add(mmd);
    }
}


void
avtPVLDFileFormat::AddBeamMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasBeam() )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::beam_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 1;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfBeamBlocks();
        md->Add(mmd);
        debug1 << "number of beam blocks: " << preader_->GetNumOfBeamBlocks() << "\n";
    }
}


void
avtPVLDFileFormat::AddSurfaceMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasSurface() )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::surface_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 2;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfSurfaceBlocks();
        md->Add(mmd);
    }
}


void
avtPVLDFileFormat::AddSphMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasSPH() )
    {
        avtMeshMetaData *mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::sph_name;
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        mmd->meshType = AVT_POINT_MESH;
        mmd->numBlocks = preader_->GetNumOfSphBlocks();
        md->Add(mmd);
    }
    else
    {
        if( add_missing_parts_  &&
                !preader_->GetMissingSphMaterial().empty() )
        {
            avtMeshMetaData *mmd = new avtMeshMetaData;
            mmd->name = PVLD_Reader::sph_name;
            mmd->spatialDimension = 3;
            mmd->topologicalDimension = 0;
            mmd->meshType = AVT_POINT_MESH;
            mmd->numBlocks = PVLD_Reader::default_number_of_partitions;
            md->Add(mmd);
        }
    }
}


void
avtPVLDFileFormat::AddTiedSetMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasTiedSet() )
    {
        const string& bn = PVLD_Reader::tiedset_name;
        avtMeshMetaData *mmd;
        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Segment" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 2;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfTiedSetBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Master" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfTiedSetBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Slave" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfTiedSetBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Bond" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 1;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfTiedSetBlocks();
        md->Add(mmd);
    }
}


void
avtPVLDFileFormat::AddContactMesh( avtDatabaseMetaData *md )
{
    if( preader_->HasContact() )
    {
        const string& bn = PVLD_Reader::contact_name;
        avtMeshMetaData *mmd;
        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Segment" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 2;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfContactBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Master" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfContactBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Slave" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfContactBlocks();
        md->Add(mmd);

        mmd = new avtMeshMetaData;
        mmd->name = PVLD_Reader::ComposeNames( bn, "Bond" );
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 1;
        mmd->meshType = AVT_UNSTRUCTURED_MESH;
        mmd->numBlocks = preader_->GetNumOfContactBlocks();
        md->Add(mmd);
    }
}




void
avtPVLDFileFormat::AddSolidVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasSolid() )
    {
        avtCentering zc = AVT_ZONECENT;
        const string& meshname = PVLD_Reader::solid_name;

        for( int i=0; i<preader_->GetNumOfSolidVariables(); i++ )
        {
            const string& name = preader_->GetSolidVariableName(i);
            const vector<int>& dims = preader_->GetSolidVariableDims(i);

            if( name == "Nodes" ) continue; // this is connectivity
            if( name == PVLD_Reader::number_of_history_name ) continue;
            if( name == PVLD_Reader::history_name ) continue;

            string varname = PVLD_Reader::ComposeNames( meshname, name );

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
            else if( dims.size()==2 )
            {
                if( dims[1]<=3 )
                {
                    avtVectorMetaData *vmd = new avtVectorMetaData;
                    vmd->name = varname;
                    vmd->meshName = meshname;
                    vmd->centering = zc;
                    vmd->hasUnits = false;
                    //vmd->varDim = (int)dims[1]; // ??? is it necessary
                    md->Add(vmd);

                    for( int d=0; d<dims[1]; d++ )
                    {
                        char buf[100];
                        sprintf(buf,"%d",d);
                        Expression *exp= new Expression;
                        exp->SetName( varname+"_"+buf );
                        exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                        exp->SetType( Expression::ScalarMeshVar );
                        exp->SetHidden( false );
                        md->AddExpression(exp);
                        delete exp;
                    }
                }
                else
                {
                    avtTensorMetaData *tensor = new avtTensorMetaData();
                    tensor->name = varname;
                    tensor->meshName = meshname;
                    tensor->centering = zc;
                    tensor->hasUnits = false;
                    //tensor->dim = 9; //(int)dims[1];
                    md->Add(tensor);

                    if( name == "Stress" )
                    {
                        string orgname = PVLD_Reader::ComposeNames( meshname, name );

                        Expression *pre= new Expression;
                        string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                        pre->SetName( prename );
                        pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                        pre->SetType( Expression::ScalarMeshVar );
                        pre->SetHidden( false );
                        md->AddExpression(pre);
                        delete pre;

                        Expression *pst= new Expression;
                        string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                        pst->SetName( pstname );
                        pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                        pst->SetType( Expression::VectorMeshVar );
                        pst->SetHidden( false );
                        md->AddExpression(pst);
                        delete pst;

                        Expression *vms= new Expression;
                        string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                        vms->SetName( vmsname );
                        vms->SetDefinition( "sqrt( 0.5*( (<" +
                                            pstname + ">[0]-<" +
                                            pstname + ">[1])^2 + (<" +
                                            pstname + ">[1]-<" +
                                            pstname + ">[2])^2 + (<" +
                                            pstname + ">[2]-<" +
                                            pstname + ">[0])^2 ))");
                        vms->SetType( Expression::ScalarMeshVar );
                        vms->SetHidden( false );
                        md->AddExpression(vms);
                        delete vms;
                    }
                }
            }
        }

        if( preader_->HasSolidHistoryVariables() )
        {
            int mx = preader_->MaxSolidHistoryVariables();
            for( int i=0; i<mx; i++ )
            {
                std::stringstream buf;
                buf << PVLD_Reader::history_name << "_" << (i+1);
                string name = buf.str();
                string varname = PVLD_Reader::ComposeNames( meshname, name );

                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }

        if( preader_->HasNode() ) AddNodeVariables( meshname, md );
    }
}



void
avtPVLDFileFormat::AddShellVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasShell() )
    {
        avtCentering zc = AVT_ZONECENT;
        const string& meshname = PVLD_Reader::shell_name;

        for( int i=0; i<preader_->GetNumOfShellVariables(); i++ )
        {
            const string& name = preader_->GetShellVariableName(i);
            const vector<int>& dims = preader_->GetShellVariableDims(i);

            if( name == "Nodes" ) continue; // this is connectivity
            if( name == PVLD_Reader::number_of_history_name ) continue;
            if( name == PVLD_Reader::history_name ) continue;

            string varname = PVLD_Reader::ComposeNames( meshname, name );

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
            else if( dims.size()==2 )
            {
                if( dims[1]<=3 )
                {
                    avtVectorMetaData *vmd = new avtVectorMetaData;
                    vmd->name = varname;
                    vmd->meshName = meshname;
                    vmd->centering = zc;
                    vmd->hasUnits = false;
                    //vmd->varDim = (int)dims[1]; // ??? is it necessary
                    md->Add(vmd);

                    for( int d=0; d<dims[1]; d++ )
                    {
                        char buf[100];
                        sprintf(buf,"%d",d);
                        Expression *exp= new Expression;
                        exp->SetName( varname+"_"+buf );
                        exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                        exp->SetType( Expression::ScalarMeshVar );
                        exp->SetHidden( false );
                        md->AddExpression(exp);
                        delete exp;
                    }
                }
                else
                {
                    avtTensorMetaData *tensor = new avtTensorMetaData();
                    tensor->name = varname;
                    tensor->meshName = meshname;
                    tensor->centering = zc;
                    tensor->hasUnits = false;
                    //tensor->dim = 9; //(int)dims[1];
                    md->Add(tensor);

                    if( name == "Stress" )
                    {
                        string orgname = PVLD_Reader::ComposeNames( meshname, name );

                        Expression *pre= new Expression;
                        string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                        pre->SetName( prename );
                        pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                        pre->SetType( Expression::ScalarMeshVar );
                        pre->SetHidden( false );
                        md->AddExpression(pre);
                        delete pre;

                        Expression *pst= new Expression;
                        string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                        pst->SetName( pstname );
                        pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                        pst->SetType( Expression::VectorMeshVar );
                        pst->SetHidden( false );
                        md->AddExpression(pst);
                        delete pst;

                        Expression *vms= new Expression;
                        string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                        vms->SetName( vmsname );
                        vms->SetDefinition( "sqrt( 0.5*( (<" +
                                            pstname + ">[0]-<" +
                                            pstname + ">[1])^2 + (<" +
                                            pstname + ">[1]-<" +
                                            pstname + ">[2])^2 + (<" +
                                            pstname + ">[2]-<" +
                                            pstname + ">[0])^2 ))");
                        vms->SetType( Expression::ScalarMeshVar );
                        vms->SetHidden( false );
                        md->AddExpression(vms);
                        delete vms;
                    }
                }
            }
        }

        if( preader_->HasShellHistoryVariables() )
        {
            int mx = preader_->MaxShellHistoryVariables();
            for( int i=0; i<mx; i++ )
            {
                std::stringstream buf;
                buf << PVLD_Reader::history_name << "_" << (i+1);
                string name = buf.str();
                string varname = PVLD_Reader::ComposeNames( meshname, name );

                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }

        if( preader_->HasNode() ) AddNodeVariables( meshname, md );
    }
}



void
avtPVLDFileFormat::AddBeamVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasBeam() )
    {
        avtCentering zc = AVT_ZONECENT;
        const string & meshname = PVLD_Reader::beam_name;

        for( int i=0; i<preader_->GetNumOfBeamVariables(); i++ )
        {
            const string& name = preader_->GetBeamVariableName(i);
            const vector<int>& dims = preader_->GetBeamVariableDims(i);

            if( name == "Nodes" ) continue; // this is connectivity
            if( name == PVLD_Reader::number_of_history_name ) continue;
            if( name == PVLD_Reader::history_name ) continue;

            string varname = PVLD_Reader::ComposeNames( meshname, name );

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
            else if( dims.size()==2 )
            {
                if( dims[1]<=3 )
                {
                    avtVectorMetaData *vmd = new avtVectorMetaData;
                    vmd->name = varname;
                    vmd->meshName = meshname;
                    vmd->centering = zc;
                    vmd->hasUnits = false;
                    //vmd->varDim = (int)dims[1]; // ??? is it necessary
                    md->Add(vmd);

                    for( int d=0; d<dims[1]; d++ )
                    {
                        char buf[100];
                        sprintf(buf,"%d",d);
                        Expression *exp= new Expression;
                        exp->SetName( varname+"_"+buf );
                        exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                        exp->SetType( Expression::ScalarMeshVar );
                        exp->SetHidden( false );
                        md->AddExpression(exp);
                        delete exp;
                    }
                }
                else
                {
                    avtTensorMetaData *tensor = new avtTensorMetaData();
                    tensor->name = varname;
                    tensor->meshName = meshname;
                    tensor->centering = zc;
                    tensor->hasUnits = false;
                    //tensor->dim = 9; //(int)dims[1];
                    md->Add(tensor);

                    if( name == "Stress" )
                    {
                        string orgname = PVLD_Reader::ComposeNames( meshname, name );

                        Expression *pre= new Expression;
                        string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                        pre->SetName( prename );
                        pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                        pre->SetType( Expression::ScalarMeshVar );
                        pre->SetHidden( false );
                        md->AddExpression(pre);
                        delete pre;

                        Expression *pst= new Expression;
                        string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                        pst->SetName( pstname );
                        pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                        pst->SetType( Expression::VectorMeshVar );
                        pst->SetHidden( false );
                        md->AddExpression(pst);
                        delete pst;

                        Expression *vms= new Expression;
                        string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                        vms->SetName( vmsname );
                        vms->SetDefinition( "sqrt( 0.5*( (<" +
                                            pstname + ">[0]-<" +
                                            pstname + ">[1])^2 + (<" +
                                            pstname + ">[1]-<" +
                                            pstname + ">[2])^2 + (<" +
                                            pstname + ">[2]-<" +
                                            pstname + ">[0])^2 ))");
                        vms->SetType( Expression::ScalarMeshVar );
                        vms->SetHidden( false );
                        md->AddExpression(vms);
                        delete vms;
                    }
                }
            }
        }

        if( preader_->HasBeamHistoryVariables() )
        {
            int mx = preader_->MaxBeamHistoryVariables();
            for( int i=0; i<mx; i++ )
            {
                std::stringstream buf;
                buf << PVLD_Reader::history_name << "_" << (i+1);
                string name = buf.str();
                string varname = PVLD_Reader::ComposeNames( meshname, name );

                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }

        if( preader_->HasNode() ) AddNodeVariables( meshname, md );
    }
}



void
avtPVLDFileFormat::AddSurfaceVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasSurface() )
    {
        avtCentering zc = AVT_ZONECENT;
        const string& meshname = PVLD_Reader::surface_name;

        for( int i=0; i<preader_->GetNumOfSurfaceVariables(); i++ )
        {
            const string& name = preader_->GetSurfaceVariableName(i);
            const vector<int>& dims = preader_->GetSurfaceVariableDims(i);

            if( name == "Nodes" ) continue; // this is connectivity

            string varname = PVLD_Reader::ComposeNames( meshname, name );

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
            else if( dims.size()==2 )
            {
                if( dims[1]<=3 )
                {
                    avtVectorMetaData *vmd = new avtVectorMetaData;
                    vmd->name = varname;
                    vmd->meshName = meshname;
                    vmd->centering = zc;
                    vmd->hasUnits = false;
                    //vmd->varDim = (int)dims[1]; // ??? is it necessary
                    md->Add(vmd);

                    for( int d=0; d<dims[1]; d++ )
                    {
                        char buf[100];
                        sprintf(buf,"%d",d);
                        Expression *exp= new Expression;
                        exp->SetName( varname+"_"+buf );
                        exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                        exp->SetType( Expression::ScalarMeshVar );
                        exp->SetHidden( false );
                        md->AddExpression(exp);
                        delete exp;
                    }
                }
                else
                {
                    avtTensorMetaData *tensor = new avtTensorMetaData();
                    tensor->name = varname;
                    tensor->meshName = meshname;
                    tensor->centering = zc;
                    tensor->hasUnits = false;
                    //tensor->dim = 9; //(int)dims[1];
                    md->Add(tensor);

                    if( name == "Stress" )
                    {
                        string orgname = PVLD_Reader::ComposeNames( meshname, name );

                        Expression *pre= new Expression;
                        string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                        pre->SetName( prename );
                        pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                        pre->SetType( Expression::ScalarMeshVar );
                        pre->SetHidden( false );
                        md->AddExpression(pre);
                        delete pre;

                        Expression *pst= new Expression;
                        string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                        pst->SetName( pstname );
                        pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                        pst->SetType( Expression::VectorMeshVar );
                        pst->SetHidden( false );
                        md->AddExpression(pst);
                        delete pst;

                        Expression *vms= new Expression;
                        string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                        vms->SetName( vmsname );
                        vms->SetDefinition( "sqrt( 0.5*( (<" +
                                            pstname + ">[0]-<" +
                                            pstname + ">[1])^2 + (<" +
                                            pstname + ">[1]-<" +
                                            pstname + ">[2])^2 + (<" +
                                            pstname + ">[2]-<" +
                                            pstname + ">[0])^2 ))");
                        vms->SetType( Expression::ScalarMeshVar );
                        vms->SetHidden( false );
                        md->AddExpression(vms);
                        delete vms;
                    }
                }
            }
        }

        if( preader_->HasNode() ) AddNodeVariables( meshname, md );
    }
}



void
avtPVLDFileFormat::AddNodeVariables( const string& meshname, avtDatabaseMetaData *md )
{
    avtCentering zc = AVT_NODECENT;

    for( int i=0; i<preader_->GetNumOfNodeVariables(); i++ )
    {
        string name = preader_->GetNodeVariableName(i);
        const vector<int>& dims = preader_->GetNodeVariableDims(i);

        //if( name == "Coordinate" ) continue; // this is connectivity
        if( name == "Index" ) name = PVLD_Reader::node_index_name;

        string varname = PVLD_Reader::ComposeNames( meshname, name );

        if( dims.size()==1 )
        {
            avtScalarMetaData *smd = new avtScalarMetaData;
            smd->name = varname;
            smd->meshName = meshname;
            smd->centering = zc;
            smd->hasUnits = false;
            md->Add(smd);
        }
        else if( dims.size()==2 )
        {
            if( dims[1]<=3 )
            {
                avtVectorMetaData *vmd = new avtVectorMetaData;
                vmd->name = varname;
                vmd->meshName = meshname;
                vmd->centering = zc;
                vmd->hasUnits = false;
                //vmd->varDim = (int)dims[1]; // ??? is it necessary
                md->Add(vmd);

                for( int d=0; d<dims[1]; d++ )
                {
                    char buf[100];
                    sprintf(buf,"%d",d);
                    Expression *exp= new Expression;
                    exp->SetName( varname+"_"+buf );
                    exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                    exp->SetType( Expression::ScalarMeshVar );
                    exp->SetHidden( false );
                    md->AddExpression(exp);
                    delete exp;
                }
            }
            else
            {
                avtTensorMetaData *tensor = new avtTensorMetaData();
                tensor->name = varname;
                tensor->meshName = meshname;
                tensor->centering = zc;
                tensor->hasUnits = false;
                //tensor->dim = 9; //(int)dims[1];
                md->Add(tensor);

                if( name == "Stress" )
                {
                    string orgname = PVLD_Reader::ComposeNames( meshname, name );

                    Expression *pre= new Expression;
                    string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                    pre->SetName( prename );
                    pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                    pre->SetType( Expression::ScalarMeshVar );
                    pre->SetHidden( false );
                    md->AddExpression(pre);
                    delete pre;

                    Expression *pst= new Expression;
                    string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                    pst->SetName( pstname );
                    pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                    pst->SetType( Expression::VectorMeshVar );
                    pst->SetHidden( false );
                    md->AddExpression(pst);
                    delete pst;

                    Expression *vms= new Expression;
                    string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                    vms->SetName( vmsname );
                    vms->SetDefinition( "sqrt( 0.5*( (<" +
                                        pstname + ">[0]-<" +
                                        pstname + ">[1])^2 + (<" +
                                        pstname + ">[1]-<" +
                                        pstname + ">[2])^2 + (<" +
                                        pstname + ">[2]-<" +
                                        pstname + ">[0])^2 ))");
                    vms->SetType( Expression::ScalarMeshVar );
                    vms->SetHidden( false );
                    md->AddExpression(vms);
                    delete vms;
                }
            }
        }
    }
}




void
avtPVLDFileFormat::AddSphVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasSPH() ||
            ( add_missing_parts_  &&
              !preader_->GetMissingSphMaterial().empty() ) )
    {
        avtCentering zc = AVT_ZONECENT;
        const string& meshname = PVLD_Reader::sph_name;

        for( int i=0; i<preader_->GetNumOfSphVariables(); i++ )
        {
            const string& name = preader_->GetSphVariableName(i);
            const vector<int>& dims = preader_->GetSphVariableDims(i);

            if( name == "Coordinate" ) continue; // this is connectivity
            if( name == PVLD_Reader::number_of_history_name ) continue;
            if( name == PVLD_Reader::history_name ) continue;

            string varname = PVLD_Reader::ComposeNames( meshname, name );

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
            else if( dims.size()==2 )
            {
                if( dims[1]<=3 )
                {
                    avtVectorMetaData *vmd = new avtVectorMetaData;
                    vmd->name = varname;
                    vmd->meshName = meshname;
                    vmd->centering = zc;
                    vmd->hasUnits = false;
                    //vmd->varDim = (int)dims[1]; // ??? is it necessary
                    md->Add(vmd);

                    for( int d=0; d<dims[1]; d++ )
                    {
                        char buf[100];
                        sprintf(buf,"%d",d);
                        Expression *exp= new Expression;
                        exp->SetName( varname+"_"+buf );
                        exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                        exp->SetType( Expression::ScalarMeshVar );
                        exp->SetHidden( false );
                        md->AddExpression(exp);
                        delete exp;
                    }
                }
                else
                {
                    avtTensorMetaData *tensor = new avtTensorMetaData();
                    tensor->name = varname;
                    tensor->meshName = meshname;
                    tensor->centering = zc;
                    tensor->hasUnits = false;
                    //tensor->dim = 9; //(int)dims[1];
                    md->Add(tensor);

                    if( name == "Stress" )
                    {
                        string orgname = PVLD_Reader::ComposeNames( meshname, name );

                        Expression *pre= new Expression;
                        string prename = PVLD_Reader::ComposeNames( meshname,"Pressure" );
                        pre->SetName( prename );
                        pre->SetDefinition( "-0.33333*trace(<" + orgname + ">)" );
                        pre->SetType( Expression::ScalarMeshVar );
                        pre->SetHidden( false );
                        md->AddExpression(pre);
                        delete pre;

                        Expression *pst= new Expression;
                        string pstname = PVLD_Reader::ComposeNames( meshname,"Principal_stress" );
                        pst->SetName( pstname );
                        pst->SetDefinition( "principal_tensor(<" + orgname + ">)" );
                        pst->SetType( Expression::VectorMeshVar );
                        pst->SetHidden( false );
                        md->AddExpression(pst);
                        delete pst;

                        Expression *vms= new Expression;
                        string vmsname = PVLD_Reader::ComposeNames( meshname,"von_Mises_Criterion" );
                        vms->SetName( vmsname );
                        vms->SetDefinition( "sqrt( 0.5*( (<" +
                                            pstname + ">[0]-<" +
                                            pstname + ">[1])^2 + (<" +
                                            pstname + ">[1]-<" +
                                            pstname + ">[2])^2 + (<" +
                                            pstname + ">[2]-<" +
                                            pstname + ">[0])^2 ))");
                        vms->SetType( Expression::ScalarMeshVar );
                        vms->SetHidden( false );
                        md->AddExpression(vms);
                        delete vms;
                    }
                }
            }
        }

        if( preader_->HasSphHistoryVariables() )
        {
            int mx = preader_->MaxSphHistoryVariables();
            for( int i=0; i<mx; i++ )
            {
                std::stringstream buf;
                buf << PVLD_Reader::history_name << "_" << (i+1);
                string name = buf.str();
                string varname = PVLD_Reader::ComposeNames( meshname, name );

                avtScalarMetaData *smd = new avtScalarMetaData;
                smd->name = varname;
                smd->meshName = meshname;
                smd->centering = zc;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }
    }
    // else {
    //   if( add_missing_parts_  &&
    //     !(preader_->GetMissingSphMaterial().empty()) ) {
    //     avtCentering zc = AVT_NODECENT;
    //     auto& meshname = PVLD_Reader::sph_name;
    // }
}




void
avtPVLDFileFormat::AddTiedSetVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasTiedSet() )
    {
        const string& basename = PVLD_Reader::tiedset_name;
        string mesh1name = PVLD_Reader::ComposeNames( basename, "Segment" );
        string mesh2name = PVLD_Reader::ComposeNames( basename, "Master" );
        string mesh3name = PVLD_Reader::ComposeNames( basename, "Slave" );
        string mesh4name = PVLD_Reader::ComposeNames( basename, "Bond" );

        for( int i=0; i<preader_->GetNumOfTiedSetVariables(); i++ )
        {
            const string& name = preader_->GetTiedSetVariableName(i);
            const vector<int>& dims = preader_->GetTiedSetVariableDims(i);

            if( name == "MstSeg"  ) continue;
            if( name == "XMaster" ) continue;
            if( name == "XSlave"  ) continue;

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd;
                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh1name,  name );
                smd->meshName = mesh1name;
                smd->centering = AVT_ZONECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh2name,  name );
                smd->meshName = mesh2name;
                smd->centering = AVT_NODECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh3name,  name );
                smd->meshName = mesh3name;
                smd->centering = AVT_NODECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh4name,  name );
                smd->meshName = mesh4name;
                smd->centering = AVT_ZONECENT;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }
    }
}


void
avtPVLDFileFormat::AddContactVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasContact() )
    {
        const string& basename = PVLD_Reader::contact_name;
        string mesh1name = PVLD_Reader::ComposeNames( basename, "Segment" );
        string mesh2name = PVLD_Reader::ComposeNames( basename, "Master" );
        string mesh3name = PVLD_Reader::ComposeNames( basename, "Slave" );
        string mesh4name = PVLD_Reader::ComposeNames( basename, "Bond" );

        for( int i=0; i<preader_->GetNumOfContactVariables(); i++ )
        {
            const string& name = preader_->GetContactVariableName(i);
            const vector<int>& dims = preader_->GetContactVariableDims(i);

            if( name == "MstSeg"  ) continue;
            if( name == "XMaster" ) continue;
            if( name == "XSlave"  ) continue;

            if( dims.size()==1 )
            {
                avtScalarMetaData *smd;
                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh1name,  name );
                smd->meshName = mesh1name;
                smd->centering = AVT_ZONECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh2name,  name );
                smd->meshName = mesh2name;
                smd->centering = AVT_NODECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh3name,  name );
                smd->meshName = mesh3name;
                smd->centering = AVT_NODECENT;
                smd->hasUnits = false;
                md->Add(smd);

                smd = new avtScalarMetaData;
                smd->name = PVLD_Reader::ComposeNames( mesh4name,  name );
                smd->meshName = mesh4name;
                smd->centering = AVT_ZONECENT;
                smd->hasUnits = false;
                md->Add(smd);
            }
        }
    }
}




void
avtPVLDFileFormat::
AddSolidMaterial( avtDatabaseMetaData* md )
{
    if( preader_->HasSolid() )
    {
        TRY
        {
            //preader_->ReadSolidMaterial();
            vector<int> mat(preader_->GetSolidMaterial());
            const vector<int>& msm = preader_->GetMissingSolidMaterial();
            for( vector<int>::const_iterator i=msm.begin(); i!=msm.end(); ++i )
                mat.push_back( *i );
            std::sort( mat.begin(), mat.end() );

            if( mat.size()>0 )
            {
                const string& meshname = PVLD_Reader::solid_name;
                avtMaterialMetaData *matmd = new avtMaterialMetaData;
                matmd->name = PVLD_Reader::ComposeNames( meshname, "mat", "_" );
                matmd->meshName = meshname;
                matmd->numMaterials = (int)mat.size();
                // std::for_each( mat.begin(), mat.end(),
                //            [&]( const int& m ) {
                //          auto& mt = preader_->GetMaterialName( m );
                //          matmd->materialNames.push_back( mt );
                //            } );
                for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                {
                    const int& m = *iter;
                    const string& mt = preader_->GetMaterialName( m );
                    matmd->materialNames.push_back( mt );
                }
                md->Add(matmd);
            }
        }
        CATCH2( std::exception, e )
        {
            string msg = e.what();
            msg += "Failed in avtPVLDFileFormat::AddSolidMaterial()\n";
            throw std::runtime_error( msg );
        }
        ENDTRY
    }
}


void
avtPVLDFileFormat::
AddShellMaterial( avtDatabaseMetaData* md )
{
    if( preader_->HasShell() )
    {
        TRY
        {
            //preader_->ReadShellMaterial();
            vector<int> mat(preader_->GetShellMaterial());
            const vector<int>& msm = preader_->GetMissingShellMaterial();
            for( vector<int>::const_iterator i=msm.begin(); i!=msm.end(); ++i )
                mat.push_back( *i );
            std::sort( mat.begin(), mat.end() );

            if( mat.size()>0 )
            {
                const string& meshname = PVLD_Reader::shell_name;
                avtMaterialMetaData *matmd = new avtMaterialMetaData;
                matmd->name = PVLD_Reader::ComposeNames( meshname, "mat", "_" );
                matmd->meshName = meshname;
                matmd->numMaterials = (int)mat.size();
                // std::for_each( mat.begin(), mat.end(),
                //            [&]( const int& m ) {
                //          auto& mt = preader_->GetMaterialName( m );
                //          matmd->materialNames.push_back( mt );
                //            } );
                for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                {
                    const int& m = *iter;
                    const string& mt = preader_->GetMaterialName( m );
                    matmd->materialNames.push_back( mt );
                }
                md->Add(matmd);
            }
        }
        CATCH2( std::exception, e )
        {
            string msg = e.what();
            msg += "Failed in avtPVLDFileFormat::AddShellMaterial()\n";
            throw std::runtime_error( msg );
        }
        ENDTRY
    }
}


void
avtPVLDFileFormat::
AddBeamMaterial( avtDatabaseMetaData* md )
{
    if( preader_->HasBeam() )
    {
        TRY
        {
            //preader_->ReadBeamMaterial();
            vector<int> mat(preader_->GetBeamMaterial());
            const vector<int>& msm = preader_->GetMissingBeamMaterial();
            for( vector<int>::const_iterator i=msm.begin(); i!=msm.end(); ++i )
                mat.push_back( *i );
            std::sort( mat.begin(), mat.end() );

            if( mat.size()>0 )
            {
                const string& meshname = PVLD_Reader::beam_name;
                avtMaterialMetaData *matmd = new avtMaterialMetaData;
                matmd->name = PVLD_Reader::ComposeNames( meshname, "mat", "_" );
                matmd->meshName = meshname;
                matmd->numMaterials = (int)mat.size();
                // std::for_each( mat.begin(), mat.end(),
                //            [&]( const int& m ) {
                //          auto& mt = preader_->GetMaterialName( m );
                //          matmd->materialNames.push_back( mt );
                //            } );
                for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                {
                    const int& m = *iter;
                    const string& mt = preader_->GetMaterialName( m );
                    matmd->materialNames.push_back( mt );
                }
                md->Add(matmd);
            }
        }
        CATCH2( std::exception, e )
        {
            string msg = e.what();
            msg += "Failed in avtPVLDFileFormat::AddBeamMaterial()\n";
            throw std::runtime_error( msg );
        }
        ENDTRY
    }
}


void
avtPVLDFileFormat::
AddSurfaceMaterial( avtDatabaseMetaData* md )
{
    if( preader_->HasSurface() )
    {
        TRY
        {
            preader_->ReadSurfaceMaterial();
            vector<int> mat = preader_->GetSurfaceMaterial();
            std::sort( mat.begin(), mat.end() );
            if( mat.size()>0 )
            {
                const string& meshname = PVLD_Reader::surface_name;
                avtMaterialMetaData *matmd = new avtMaterialMetaData;
                matmd->name = PVLD_Reader::ComposeNames( meshname, "mat", "_" );
                matmd->meshName = meshname;
                matmd->numMaterials = (int)mat.size();
                // std::for_each( mat.begin(), mat.end(),
                //            [&]( const int& m ) {
                //          auto& mt = preader_->GetMaterialName( m );
                //          matmd->materialNames.push_back( mt );
                //            } );
                for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                {
                    const int& m = *iter;
                    const string& mt = preader_->GetMaterialName( m );
                    matmd->materialNames.push_back( mt );
                }
                md->Add(matmd);
            }
        }
        CATCH2( std::exception, e )
        {
            string msg = e.what();
            msg += "Failed in avtPVLDFileFormat::AddSurfaceMaterial()\n";
            throw std::runtime_error( msg );
        }
        ENDTRY
    }
}


void
avtPVLDFileFormat::
AddSphMaterial( avtDatabaseMetaData* md )
{
    if( preader_->HasSPH() ||
            ( add_missing_parts_  &&
              !preader_->GetMissingSphMaterial().empty() ) )
    {
        TRY
        {
            //preader_->ReadSphMaterial();
            vector<int> mat(preader_->GetSphMaterial());
            const vector<int> &msm = preader_->GetMissingSphMaterial();
            for( vector<int>::const_iterator i=msm.begin(); i!=msm.end(); ++i )
                mat.push_back( *i );
            std::sort( mat.begin(), mat.end() );

            if( mat.size()>0 )
            {
                const string& meshname = PVLD_Reader::sph_name;
                avtMaterialMetaData *matmd = new avtMaterialMetaData;
                matmd->name = PVLD_Reader::ComposeNames( meshname, "mat", "_" );
                matmd->meshName = meshname;
                matmd->numMaterials = (int)mat.size();
                // std::for_each( mat.begin(), mat.end(),
                //            [&]( const int& m ) {
                //          auto& mt = preader_->GetMaterialName( m );
                //          matmd->materialNames.push_back( mt );
                //            } );
                for( vector<int>::const_iterator iter=mat.begin(); iter!=mat.end(); iter++ )
                {
                    const int& m = *iter;
                    const string& mt = preader_->GetMaterialName( m );
                    matmd->materialNames.push_back( mt );
                }
                md->Add(matmd);
            }
        }
        CATCH2( std::exception, e )
        {
            string msg = e.what();
            msg += "Failed in avtPVLDFileFormat::AddSphMaterial()\n";
            throw std::runtime_error( msg );
        }
        ENDTRY
    }
}

// ****************************************************************************
// Class: vtkFloatArrayFromFloatVector
//
// Purpose:
//   Steals a vector<float>'s data and exposes it as a vtkPoints object. The
//   vector data gets freed when the vtkPoints object gets freed.
//
// Notes:
//
// Programmer: Brad Whitlock
// Creation:   Mon Jun  2 17:32:25 PDT 2014
//
// Modifications:
//
// ****************************************************************************

class vtkFloatArrayFromFloatVector : public vtkCommand
{
public:
    static vtkFloatArrayFromFloatVector *New()
    {
        return new vtkFloatArrayFromFloatVector;
    }

    // Install the observer on the DeleteEvent.
    vtkFloatArray *MakeArray(int nComponents, vector<float> &obj)
    {
        // Steal the data from the input vector by swapping with a new vector
        // that we control.
        this->vec = new vector<float>();
        this->vec->swap(obj);

        // Wrap the stolen data as a vtkPoints. Indicate that the vtkPoints
        // object is not allowed to delete the vector storage.
        vtkFloatArray *arr = vtkFloatArray::New();
        arr->SetNumberOfComponents(nComponents);
        arr->SetArray(&this->vec->operator[](0), this->vec->size(), 1);

        // Ensure that this object will be called in response to the
        // vtkPoints object getting deleted.
        this->Id = arr->AddObserver(vtkCommand::DeleteEvent, this);

        return arr;
    }

    // This method gets invoked when the vtkPoints object gets deleted.
    virtual void Execute(vtkObject *object, unsigned long, void *)
    {
        // Delete the vector data.
        delete this->vec;
        this->vec = NULL;

        // Remove this object as an observer on the vtkPoints.
        object->RemoveObserver(this->Id);

        // Delete this observer object.
        this->UnRegister();
    }

protected:
    vtkFloatArrayFromFloatVector() : vec(NULL), Id(0) {}
    virtual ~vtkFloatArrayFromFloatVector() {}

private:
    vector<float> *vec;
    unsigned long  Id;

    void operator=(const vtkFloatArrayFromFloatVector &); // not implemented
    vtkFloatArrayFromFloatVector(const vtkFloatArrayFromFloatVector &); // not implemented
};

// ****************************************************************************
// Method: avtPVLDFileFormat::MakePoints
//
// Purpose:
//   Turns the input coordinates into vtkPoints.
//
// Arguments:
//   crd : The vector containing the input coordinates.
//
// Returns:    A new vtkPoints object that steals the vector's data.
//
// Note:       The active implementation creates a new vector, swaps storage
//             with the input vector, and wraps that new vector's array as
//             vtkPoints. The storage is deleted later when the vtkPoints object
//             gets deleted.
//
// Programmer: Brad Whitlock
// Creation:   Mon Jun  2 17:30:22 PDT 2014
//
// Modifications:
//
// ****************************************************************************

vtkPoints *
avtPVLDFileFormat::MakePoints(vector<float> &crd)
{
    vtkPoints *pts = vtkPoints::New();

    if(enablePVLD)
    {
        // Note: we create this object to steal the array from the vector. The vector
        // storage gets deleted when the vtkPoints is deleted -- and so does the
        // vtkFloatArrayFromFloatVector instance.
        vtkFloatArray *arr = vtkFloatArrayFromFloatVector::New()->MakeArray(3, crd);
        pts->SetData(arr);
        arr->Delete();
    }
    else
    {
        int nvrt = (int)crd.size() / 3;
        pts->SetNumberOfPoints(nvrt);
        float *pt = (float *) pts->GetVoidPointer(0);
        for( size_t i=0; i<crd.size(); i++ )
            pt[i] = crd[i];
    }
    return pts;
}

// ****************************************************************************
// Method: avtPVLDFileFormat::GenerateSolidMesh
//
// Purpose:
//   Create an unstructured mesh for the solid.
//
// Arguments:
//   crd : The coordinates.
//   ele : The elements.
//
// Returns:    A vtkUnstructuredGrid containing the solid mesh.
//
// Note:       crd's data gets stolen so we don't have to create a copy.
//
// Programmer: hao pan
// Creation:   Mon Jun  2 17:29:32 PDT 2014
//
// Modifications:
//    Brad Whitlock, Mon Jun  2 17:29:50 PDT 2014
//    Made it use MakePoints, which steals the storage from the vector and
//    wraps it as vtkPoints.
//
// ****************************************************************************

vtkUnstructuredGrid*
avtPVLDFileFormat::GenerateSolidMesh( vector<float>& crd, const vector<int>& elm )
{
    int nvrt = (int)crd.size()/3; (void) nvrt;
    int nele = (int)elm.size()/8;

    vtkPoints *pts = MakePoints(crd);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate(nele);
    pts->Delete();

    vtkIdType verts[8];
    const int* cnt=&elm[0];
    for( int i=0; i<nele; i++ )
    {
        for( int j=0; j<8; j++ )
            verts[j] = *cnt++;

        if( verts[3] == verts[4] )   // tetra
            ugrid->InsertNextCell( VTK_TETRA, 4, verts );
        else if( verts[4] == verts[5]  && verts[6] == verts[7]  )   // prism
        {
            vtkIdType prism[6];
            prism[0]=verts[0];
            prism[1]=verts[1];
            prism[2]=verts[4];
            prism[3]=verts[3];
            prism[4]=verts[2];
            prism[5]=verts[6];
            ugrid->InsertNextCell( VTK_WEDGE, 6, prism );
        }
        else
            ugrid->InsertNextCell( VTK_HEXAHEDRON, 8, verts );
    }
    return ugrid;
}



vtkUnstructuredGrid*
avtPVLDFileFormat::GenerateShellMesh( vector<float>& crd, const vector<int>& elm )
{
    int nele = (int)elm.size()/4;

    vtkPoints *pts = MakePoints(crd);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate(nele);
    pts->Delete();

    vtkIdType verts[4];
    const int* cnt=&elm[0];
    for( int i=0; i<nele; i++ )
    {
        for( int j=0; j<4; j++ )
            verts[j] = *cnt++;
        ugrid->InsertNextCell( VTK_QUAD, 4, verts );
    }
    return ugrid;
}



vtkUnstructuredGrid*
avtPVLDFileFormat::GenerateBeamMesh( vector<float>& crd, const vector<int>& elm )
{
    int nele = (int)elm.size()/2;

    vtkPoints *pts = MakePoints(crd);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate(nele);
    pts->Delete();

    vtkIdType verts[2];
    const int* cnt=&elm[0];
    for( int i=0; i<nele; i++ )
    {
        for( int j=0; j<2; j++ )
            verts[j] = *cnt++;
        ugrid->InsertNextCell( VTK_LINE, 2, verts );
    }
    return ugrid;
}




vtkUnstructuredGrid*
avtPVLDFileFormat::GenerateSphMesh( vector<float>& crd, const vector<int>& elm )
{
    int nele = (int)elm.size();

    vtkPoints *pts = MakePoints(crd);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate(nele);
    pts->Delete();

    vtkIdType verts[1];
    const int* cnt=&elm[0];
    for( int i=0; i<nele; i++ )
    {
        verts[0] = *cnt++;
        ugrid->InsertNextCell( VTK_VERTEX, 1, verts );
    }
    return ugrid;
}




vtkUnstructuredGrid*
avtPVLDFileFormat::GeneratePointMesh( vector<float>& crd )
{
    int nvrt = (int)crd.size()/3;

    vtkPoints *pts = MakePoints(crd);
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate( nvrt );
    pts->Delete();

    vtkIdType verts[1];
    for( int i=0; i<nvrt; i++ )
    {
        verts[0] = i;
        ugrid->InsertNextCell( VTK_VERTEX, 1, verts );
    }
    return ugrid;
}


vtkUnstructuredGrid*
avtPVLDFileFormat::GenerateBondMesh( const vector<float>& crd1, const vector<float>& crd2 )
{
    int nvrt = (int)crd1.size()/3;

    vtkPoints *pts = vtkPoints::New();
    pts->SetNumberOfPoints(2*nvrt);

    float *pt = (float *) pts->GetVoidPointer(0);
    for( int i=0; i<nvrt; i++ )
    {
        pt[6*i  ] = crd1[3*i  ];
        pt[6*i+1] = crd1[3*i+1];
        pt[6*i+2] = crd1[3*i+2];
        pt[6*i+3] = crd2[3*i  ];
        pt[6*i+4] = crd2[3*i+1];
        pt[6*i+5] = crd2[3*i+2];
    }

    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->SetPoints( pts );
    ugrid->Allocate( nvrt );
    pts->Delete();

    vtkIdType verts[2];
    for( int i=0; i<nvrt; i++ )
    {
        verts[0] = 2*i;
        verts[1] = 2*i+1;
        ugrid->InsertNextCell( VTK_LINE, 2, verts );
    }
    return ugrid;
}

// ****************************************************************************
//
// Modifications:
//   Brad Whitlock, Mon Jun  2 17:57:55 PDT 2014
//   I changed the code to steal the data from the vector and return it as
//   a vtkFloatArray so we can avoid copies.
//
// ****************************************************************************

vtkFloatArray*
avtPVLDFileFormat::GenerateVariable( const vector<int>& dims, vector<float>& dat )
{
    int ncmp = 1;
    for( size_t i=1; i<dims.size(); i++ )
        ncmp *= dims[i];
    int msz  = dims[0];

    vtkFloatArray *var = NULL;
    if(enablePVLD)
    {
        var = vtkFloatArrayFromFloatVector::New()->MakeArray(ncmp, dat);
    }
    else
    {
        var = vtkFloatArray::New();
        var->SetNumberOfComponents(ncmp);
        var->SetNumberOfTuples( msz );

        float *des = (float *)var->GetVoidPointer(0);
        const float *src = &dat[0];
        for( int i=0; i<msz*ncmp; i++ )
            *des++ = *src++;
    }

    return var;
}


vtkFloatArray*
avtPVLDFileFormat::GenerateDyna3dStressTensor( const vector<int>& dims, const vector<float>& dat )
{
    int msz = dims[0];

    vtkFloatArray* var = vtkFloatArray::New();
    var->SetNumberOfComponents(9);
    var->SetNumberOfTuples( msz );

    float *des = (float *)var->GetVoidPointer(0);
    const float *src = &dat[0];
    for( int i=0; i<msz; i++ )
    {
        PVLD_Reader::ConvDyna3dStreeTensor( src, des );
        src += 6;
        des += 9;
    }
    return var;
}



vtkFloatArray*
avtPVLDFileFormat::
GenerateSphStressTensor( const vector<int>& dims, const vector<float>& dat )
{
    int msz = dims[0];

    vtkFloatArray* var = vtkFloatArray::New();
    var->SetNumberOfComponents(9);
    var->SetNumberOfTuples( msz );

    float *des = (float *)var->GetVoidPointer(0);
    const float *src = &dat[0];
    for( int i=0; i<msz; i++ )
    {
        PVLD_Reader::ConvSphStreeTensor( src, des );
        src += 6;
        des += 9;
    }
    return var;
}


// vtkFloatArray*
// avtPVLDFileFormat::GenerateTiedSetVariable( const vector<int>& dims, const vector<float>& dat )
// {
//   int ncmp = 1;
//   for( int i=1; i<dims.size(); i++ )
//     ncmp *= dims[i];
//   int msz  = dims[0];

//   vtkFloatArray* var = vtkFloatArray::New();
//   var->SetNumberOfComponents( ncmp );
//   var->SetNumberOfTuples( msz );

//   float *des = (float *)var->GetVoidPointer(0);
//   const float *src = &dat[0];
//   for( int i=0; i<msz*ncmp; i++ ) {
//     *des++ = *src++;
//   }
//   return var;
// }




void
avtPVLDFileFormat::
AddSolidIndexVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasSolid() )
    {
        const PVLD_Reader::IndexVariables& iv = preader_->GetSolidIndexVariables();
        if( iv.names_.size()>0 )
        {
            const string& meshname = PVLD_Reader::solid_name;
            for( size_t i=0; i<iv.names_.size(); i++ )
            {
                const string& idxname = iv.names_[i];
                const vector<string>& dsname  = iv.dsnames_[i];
                const vector<vector<int> >& dsdims  = iv.dsdims_[i];
                AddIndexVariables( meshname, idxname, dsname, dsdims, md );
            }
        }
    }
}


void
avtPVLDFileFormat::
AddShellIndexVariables( avtDatabaseMetaData *md )
{
    if( preader_->HasShell() )
    {
        const PVLD_Reader::IndexVariables& iv = preader_->GetShellIndexVariables();
        if( iv.names_.size()>0 )
        {
            const string& meshname = PVLD_Reader::shell_name;
            for( size_t i=0; i<iv.names_.size(); i++ )
            {
                const string& idxname = iv.names_[i];
                const vector<string>& dsname  = iv.dsnames_[i];
                const vector<vector<int> >& dsdims  = iv.dsdims_[i];
                AddIndexVariables( meshname, idxname, dsname, dsdims, md );
            }
        }
    }
}





void
avtPVLDFileFormat::
AddIndexVariables( const string& meshname,
                   const string& idxname,
                   const vector<string>& dsname,
                   const vector<vector<int> >& dsdims,
                   avtDatabaseMetaData *md )
{
    string basename = PVLD_Reader::ComposeNames( meshname, idxname );

    for( size_t i=0; i<dsname.size(); i++ )
    {
        const string& name = dsname[i];
        const vector<int>& dims = dsdims[i];

        string varname = PVLD_Reader::ComposeNames( basename, name );

        if( dims.size()==1 )
        {
            avtScalarMetaData *smd = new avtScalarMetaData;
            smd->name = varname;
            smd->meshName = meshname;
            smd->centering = AVT_ZONECENT;
            smd->hasUnits = false;
            md->Add(smd);
        }
        else if( dims.size()==2 )
        {
            if( dims[1]<=3 )
            {
                avtVectorMetaData *vmd = new avtVectorMetaData;
                vmd->name = varname;
                vmd->meshName = meshname;
                vmd->centering = AVT_ZONECENT;
                vmd->hasUnits = false;
                md->Add(vmd);

                for( int d=0; d<dims[1]; d++ )
                {
                    char buf[100];
                    sprintf(buf,"%d",d);
                    Expression *exp= new Expression;
                    exp->SetName( varname+"_"+buf );
                    exp->SetDefinition( "array_decompose(<"+varname+">,"+buf+")");
                    exp->SetType( Expression::ScalarMeshVar );
                    exp->SetHidden( false );
                    md->AddExpression(exp);
                    delete exp;
                }
            }
        }
    }
}

